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CTN \ The boundary integral method for calculating the stationary states of a quantum particle in 

nano-devices and quantum billiards is presented in detail at an elementary level. According to the 

method, wave functions inside the domain of the device or billiard are expressed in terms of line 

^ , integrals of the wave function and its normal derivative along the domain's boundary; the respective 

D ■ energy eigenvalues are obtained as the roots of Fredholm determinants. Numerical implementations 

y^^ ' of the method are described and applied to determine the energy level statistics of billiards with 

Qs , circular and stadium shapes and demonstrate the quantum mechanical characteristics of chaotic 

motion. The treatment of other examples as well as the advantages and limitations of the boundary 

integral method are discussed. 
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■^ I. INTRODUCTION 

o , 

^ \ Recent advances in nanotechnology, based on advanced crystal growth and lithogmphic techniques, have opened an 

O ' avenue to fabricate very small and clean electronic devices, known as nano-deviceaa. The charge carriers (electrons) 

'c/3 I in such devices, through gate voltages, are confined to one or two spatial dimensions. At very low temperatures, 

J>-». the spatial extent of the systems along the direction of confinement is comparable to the Fermi wavelength of the 

tG ' electrons. Quantum dots and quantum wires are examples of quasi zero- and one-dimensional nano-devices in which 

, ^~ <i: confinement of the electrons occur along all three and along two spatial directions, respectively, while in the inversion 

layer of narrow-gap semiconductor heterostructures the electrons are confined along the_direction perpendicular to 

y—i ', the layer. Quantum dots are relevant in the study of-|thc Coulomb blockade phenomenaa, while quantum wires are 

^ ' experimental realizations of so-called Luttinger liquidM. 

Cn I The motion of the electrons in a clean two-dimensional nano-device is ballistic, i.e., the electrons are scattered 

^^ . mainly by the device boundaries and not by impurities. The device boundaries, due to high precision lithography, 
I~i ' may have arbitrary shapes and are very sharp, i.e., the electrical potential changes abruptly on atomic scales. As 
f^ , a result, the behavior of such two-dimensional nano-devices, which exhibit quantum confinement in one direction 
r~^ ■ and free motion of the electrons in a finite two-dimensional domain of sub-micron size, is governed by single-electron 
Q^ , (particle) physics, and can be described theoretically by solving the corresponding Schrodingerjvave equation. Such 
^^ nano-devices can be considered as quantum mechanical analogue of classical billiard systemajl in which point like 
O particles bounce inside a two-dimensional (2D) region V delimited by the contour F. An idealized quantum billiard 
' j/2 • confines a quantum particle inside a 2D infinite potential well; the shape of the infinite well being determined by F. 
J>-j' Quantum billiards represent models of nano-devices which play an important role in modern semiconductor 

1-^ industryGI. The experimental study, via STM techniques, of quantum billiards provides a new testing-ground for 

t^ the predictions of quantum mechanicstJ. The study of quantum billiards allows one to investigate also the quantum 
J> , signatures of classical chaos. It is known that non-integrable classical systems are chaotic, i.e., the phase space tra- 
jectory of the system exhibits exponential sensitivity to the initial conditions. In the case of billiards, the chaotic 
. _ behavior is caused by the irregularities of the boundary and not by the complexity of the interaction in the system 
jrt ■ (e.g., scattering of the particle from randomly distributed impurities). Since the concept of "phase space trajectory" 
■ ■ ' loses its meaning in quantum mechanics, one can naturally ask oneself what is the quantum mechanical analogue of 
(classical) chaos, or more precisely, is there any detectable difference between the behavior of a quantum system with 
chaotic- and non-chaotic classical limit, respectively. 

The answer to these qi:gstions should be sought in the characteristics of the fiuctuations of the energy levels of the 
quantum billiard systemsEra. Thus, in order to study the physical properties of quantum billiards one needs to find 
first the corresponding energy spectrum by solving the time independent Schrodinger equation 



HMr) 



.|;V^ + V-M 



«„(r) = E„Mr) , (1.1) 



where H is the Hamiltonian of the system, V(r ) is the potential energy, and ■)/'„ is the eigenfunction corresponding 



to the energy eigenvalue En. In general, in (LI) the potential V{r) does not contain the term corresponding to the 
infinite potential well; the effect of the later is reflected by the "hard- wall" (i.e., Dirichlet) boundary conditions at the 
billiard boundary. The spectrum is discrete and the distribution of the energy levels En is determined by the form of 
the po tent ial and by the boundary conditions. 



Eq.(Ll) can be solved analytically only for very few special cases, when the system is integrable, i.e., when there 
exists, besides the energy, a second conserved physical quantity. Such examples, like a quantum particle in a rect- 
angular or circularJnfinite potential well, are discussed in most of the quantum mechanics textbookaj and in some 
recent publicationscl, as well. However, for a generic quantum billiard the energy spectrum can be determined only 
numerically, and the description of such numerical methods lacks in all widely used quantum mechanics textbooks. 

The purpose of the present article is to fill this gap by providing the reader with a self-contained and practical intro- 
duction to a powerful numerical method, known as the Boundary Integral Method (BIM), for calculating the energy 
levels of a 2D quantum system, e.g., a quantum billiard. While the BIM, sometimes also referred to as the Bo-^ 
Element Method (BEM), has been extensively used for many years for solving different engiaeeijJag-|Ptt)blems 
its application for calculating energy spectra of quantum billiards has emerged only recentlytj'E3iiJll3'E3. 

Before we embark on our presentation of the BIM for calculating energy spectra of 2D quantum systems, let us 
first mention a few other frequently used numerical methods in the same context. 

Essentially all numerical methods devised to solve the single particle Schrodinger equation (1^) can be classified in 
two groups. The methods belonging to the first group assume that one readily knows a complete set of orthonormal 
functions {0m('')} which obey the desired boundary conditions along the billiard boundary. By expanding the 
unknown energy eigenfunctions 

i'nir) = ^c„„(/)„(r), (1.2) 

m 

the Schrodinger equation is converted into the familiar system of homogeneous linear equations for the coefficients 

/ ^ \-ti71m ^nVnm) Cnm — t( , (^i.oj 

m 

Here 6nm is the Kronecker-delta (equal to 1 for n = m and zero otherwise) , and the matrix elements of the Hamiltonian 
are 

'drrn{r)H<j)rn{r). (1.4) 



Equation (O) admits non-trivial solutions (energy eigenstates or stationary states) only for those values of En (the 



energy eigenvalues) which satisfy the condition 

det\Hnm - EnSnml = 0. (1.5) 

This condition can be employed to determine the i?„'s. 

When the billiard boundary is irregular, in general, it is impossible to find analytical expressions for the functions 
4>m{r) and, therefore, the method as described fails. However, in this case one can overcome the previously mentioned 
difficulty by either performing a coordinate transformation which renders the boundary highly regular, or by extending 
the system, fitting the billiard inside a rectangle or circle along which the Dirichlet boundary conditions apply. Now a 
complete set of orthonormal functions can be easily found, but the price one pays in both cases is that the correMionding 
Hamiltonian becomes more complicated: in the fopner case the simple form of the kinetic energy is alteredEJ while 
in the latter case the potential energy is modifiedtj, i.e., V{r) = inside V and V{r) = oo (in practice a suitably 
chosen large value) outside V. 



The second class of numerical methods intended to calcula te b illiard spectra regard Eq.(LT) as a partial differential 



equation for which the general solution is formally given by (L2) for some conveniently chosen basis functions (j)mir). 
The energy eigenfunctions and eigenvalues are determined by requiring the general solution ( |l.2| ) to obey the Dirichlet 
boundary conditions along dV. 0|Lcourse, the boundary conditions can be met only for particular values of the energy, 
i.e., the energy eigenvalues. HelleiEEl used this method choosing as the basis functions plane waves, rfiile a more general 
and systematic implementation of this method in plane polar coordinates is described by SchmitEll. 

The BIM is an efficient alternative to the above mentioned two classes of methods for solving numerically the 
Schrodinger equation. We shall consider its application only for two-dimensional systems. The BIM will allow us to 



study the quantum analogue of classical chaotic systems and reveal that chaotic behavior is reflected in the spacing 
of the energy eigenvalues E„. For this purpose, the BIM is formulated in Sec. |l| and is applied, in Sec. [II, to the 
spectra of circular, stadium and generalized stadium billiards. In Sec. IV we discuss further examples to which the 
BIM can be applied. In Sec. |V| we present concluding remarks. 



II. THE BOUNDARY INTEGRAL METHOD 



Consider a quantum particle of mass m moving in a finite, simply connected region V, experiencing the potential 
V{r) and being governed by the Hamiltonian 



2m 



V{r) 



(2.1) 



The energy spectrum of the particle can be determined from the time- independent Schrodinger equation (1.1) together 



with the boundary conditions for the wave functions 4'n{f) specified on a closed curve T = dT> which delimits the 
region V. 



The Schrodinger equation (1.1) is an implicit equation for £"„ and ^„(r). This differential equation can be replaced 
by an implicit integral equation which can also serve to determine En and '0„(r). For this purpose, one introduces 
the Green's function, G{r, r'; E) of the operator E — H, defined as the solution of 



[E-H{r)]G{r,r';E) = 6{r ^ r') 



(2.2) 



where d{r — r') is the two-dimensional (5-function, _E is a complex variable , and r, r' are arbitrary points in V. 
Multiplying Eq.(l.l) by G{r,r';E), Eq.(p.2[) by ipnir), and adding the resulting equations yield 



V'„(r)<5(r - r') + {E„ - E) VA.(r)G(r, r'; E) = G(r, r'; E)H^n{r) ~ Vn(r)-ffG(r, r'; E) . (2.3) 

We consider now Eq.( |2.3| ) for E = _E„. In this case the second term on the LHS vanishes, provided that G is finite 
(i.e., has no poles) at En- A necessary (bu t no t sufficient) condit ion is that G does not obey the same boundary 
conditions as ■)/'„. Inserting the Hamiltonian (2.1) in the RHS of Eq.( [2Ti^ ) eliminates the terms containing the potential 
energy V(r) and Eq.(p.3[) becomes 



V'„(r)(5(r-r') 



2m 



[V'„(r)V2G (r, r'; En) ~ G (r, r'; En) VVn(r-)] 



(2.4) 



Recalling the identity uV'^v — V(uVw) — VuVw, valid for any differcntiablc functions u{r) and v{r), the RHS of the 
above equation can be written as a divergence 



V'„(r)5(r- 



—V • [iPn{r)^G (r, r'; En) - G (r, r'; En) V^„(r)] 
Zra 



(2.5) 



Integration with respect to r over the domain T> yields, on the LHS, 4'n{i'') since r' € V; applying Green's formula 
the integral on the RHS can be expressed as a line integral along F = dV and Eq.(2.[) becomes 



Mr') 



ds{r) [Mr)d,G (r, r'; En) - G (r, r'; E^) dAn {r)\ 



(2.6) 



2m jy 

Here ds(r) is the infinitesimal arc length along F at r e F, and the normal derivative dy is defined through 

a, = v(r) ■ Vr , (2.7) 

with vir) representing the exterior normal unit vector to F at r S F. This is the desired integral equation whic h, fo r 
nano-devices and quantum billiards, provides a simpler avenue to En and 4'ni'i') than the Schrodinger equation ( p, . l| ) . 
Note that Eq. ( |2.6| ) does not exhibit an explicit dependence on the potential function V{r); the effect of the latter is 
incorporated entirely in the Green's function G{r,r';E). 

The eigenvalues En can be obtained by noting that existence of solutions 4'n{i') implies conditions of the type ( |1.5|) . 
We will adopt a similar strategy for Eq. ( |2.6| ) and consider for this purpose the limit r' S F. In this limit Eq.(p.6|) 
becomes an implicit equation for ipnir) confined solely to the boundary F such that a condition like (1.5) can be 
postulated and exploited to determine -E„. 



The limit r' € F in (1.5) is not trivial since both the Green's function and its normal derivative are singular at 
r — r' . However, these singu larit ies are integrable in the sense of Cauchy's principal value. To demonstrate this we 
carry out the integration in (1.5) along a slightly altered contour T^ which avoids the singularity and then let the 
altered contour approach F continuously. For this purpose we define F^ = F^ U Cg, where F^ coincides with F, except 
for a portion of arc-length 2e centered about r'; C^ is a circular arc with center at r' and radius e as shown in Fig. ^, 
where r' lies inside the region delimited by F^. We consider then the integral in (L5) for lim^^Q+T^ = F. 

For Fe the integral has two contributions corresponding to F^ and C^ ■ The integration along Fg in the e ^ 0^ limit 



e(r') 




Fig. 1. Geometry of the boundary Fg in the vicinity of the point A (position vector r') where the Green's function is singular, 
is, by definition, Cauchy's principal value integral along the original contour F. We denote the integral as 



e-»0 



lim / ds{r) 



' ® ds{r) . 



(2.8) 



The contribution due to the integral along C^ depends on the type of singularity of the Green's function at r = r'. 
The integral can be calculated as shown in Appendix B. One obtains 



2m 



1 



lim^^OTT- / ds{r)[tljn{r)dt,G{r,r'-En)~G{r,r';En)di,i:n{r)] = ■;r'0n(r') 



(2.9) 



In the derivation of this formula we have implicitly assumed that there is a unique tangent to F at r', i.e., that the 
angle 8{r') in Fig. ^lis equal to n. Otherwise, according to Eq.(B4) in Appendix B, the RHS of (2.£) must be replaced 



by (6'(r')/27r) ■(/'„(r^, where 0{r') is the exterior angle between the two tangents to F at r'. 
Altogether, one obtains for ipn{r') , r' e F the integral equation 



Mr') 



h 



V i ds{r) [Mr)d,G (r, r'; E,,) ~ G (r, r'; E^) d,Mr)] 
m 7r 



(2.10) 



where one still needs to specify the boundary condition on F which involves ipn and/or its normal derivative d^ijj„ 
The boundary condition is expressed as a linear functional relation 



T[Mr),d,Mr)] = 



r e F 



(2.11) 



The actual form of the functional JF depends on the physi cal p roble m at hand, but not on the contour F. The 
energy eigenvalues En are determined by requiring that Eqs. (|2.10|) and ( p. 11 ) admit nontrivial solutions for ipn- This 
condition leads us to an equation involving functional (Fredholm) determinants of the type (1.5) which need to be 
solved by numerical means. Once E„ and the corresponding tpn and d^tpn on F are determined, the eigenfunction 
inside the domain V can be calculated using Eq.(2.6). 

Below we will demonstrate the application of the method outlined which is referred to as the Boundary Integral 
Method (BIM). The met hod is practical whenever (i) a Green's function G is available analytically and (ii) the 
boundary condition (2.11) is fairly simple; the method applies to F of arbitrary shape. 



III. BILLIARD SPECTRA VIA BIM 



Inside a billiard a particle moves freely, i.e., y = in (2.1). The Green's function defined through (2.2) is well 
known in this case and is given by 



G{r,r';E) = -i^i/W (fc |r - r'|) , 



(3.1) 



as shown in Appendix A. Here k = \/2rriE/h is the so-called wave vector; the index n is dropped since we focus in the 
following on a single eigenstate. We will also use the notation G (r, r'; k) for the Green's function. Si nce t he particle 
is confined to the billiard, its wave function ip = ^n must vanish along F and the boundary condition ( 2.11 ) takes the 
form 



V'(t) ~ , dij'tpi''') — arbitrary 



Vr e r 



(3.2) 



Inserting (3.2) in the the integral equation (2.1C) leads to 



V (bds{r)G{r,r';E)d^^/j{r) = 



(3.3) 



This integral equation admits non-trivial solutions only if the corresponding Fredholm (functional) determinant van- 
ishes, i.e., for 

det [G(r,r';£;)] = 0, (3.4) 

a condition which allows one to determine the energies En- 

Even though the analytical expression of the Green's function G is known, the Fredholm determinant (3.4) is 
difficult to evaluate for arbitrary billiard boundaries F. Below we describe more practical schemes for solving the 
integral equation (^.3]). 



A. Methods for Solving the BIE 



There are basically three different methods for solving the BIE (2.10) for the billiard problem. Before presenting 
these methods, let us first parameterize the bilhard boundary F through the arc length s e [0,£], where £ is the 
length of the billiard boundary F. Thus, the position of each point r G F is uniquely determined by s through the 
function r = r{s). It is convenient to introduce the notation 



The BIE (3^) can be recast then as 



u{s) = Unir{s)) = d^Tpnir) 



[ dsG{s,s';k)u{s) = 0, 
Jo 



(3.5) 



(3.6) 



where, for brevity, we have dropped the index n which labels the eigenstates. 

Method I. The most obvious (but not necessarily the most convenient) method of solution relies on the obser- 
vation that both the wave function and its normal derivative (i.e., u{s)) are single-valued functions and, therefore, 
m(s) must be a periodic function of s with period C. Hence, (3.6) can be expressed as a Fourier series 



u{s) 



X! Ujexp(iKjs) , 



(3.7a) 



j = -oo 



where 



K = ^, 



j=0,±l,±2,. 



(3.7b) 



and where Uj is the Fourier transform of u(s) 



1 /-^ 

— / dsu{s)exp{—iKjs) 
•-■ Jo 



(3.7c) 



By taking the Fourier transform of Eq.(3.(:) with respect to s' and using (3.7a) one obtains the system of hnear 
equations 



Yl ^y(^)"j 



j = -oo 



where 



Aij{k) = -T / ds ds'G{s,s';k)exp[~i{Kis' - Kj 
•-' Jo Jo 



(3.8) 



(3.9) 



Here the information about the biUiard boundary F is contained in the s- and s'-dependence of the Green's function 
through r = r{s) and r' — r{s'). The energy eigenvalues, expressed through k, are the solutions of the equation 



det[Ay(fc)] = 



(3.10) 



which must hold in order to render (3.8) so lvabl e. 

For an arbitrary F one cannot solve Eq.( p.lC| ) exactly. However, approximate energy eigenvalues can be obtained 



by truncating the infinite system of linear equations (3.8) at some suitably chosen wave vector K ^- Th e truncation 



implies that the Fourier components of u{s) which correspond to \Kj \ > Kc are set equal to zero in ( |3.7a ). In this case 
the relevant part of the matrix Aij becomes finite and the corresponding determinant can be calculated numerically. 
The drawback of the truncation is that the resulting energy eigenvalues expressed through k are accurate only as long 
as fc ^ Kc- If one seeks to describe energy levels with larger fc-values one needs to increase Kc which, however, leads 
to an increased computational effort, the latter increasing rapidly with the dimension of the matrix Aij. 

The calculation of the matrix elements Aij as double integrals (with an integrand which is singular at s = s') is 
computationally cumbersome and, as a result, the present method is impractical, except for the case of a circular 
billiard. In this case Aij (k) is a diagonal matrix and it s ele ments can be expressed in terms of products of Bessel and 
Hankel functions as shown in Appendix C. Equation ( 3.1C| ) reads then 



det [A,j{k)] 



ex 



n 



Mk)Hl'\k) 



0. 



(3.11) 



The Hankel functions H- (k) have no real roots and, hence, the energy eigenvalues for a circular billiard with unit 
radius are given by the zeros of the integer order Bessel functions 



Jj(fc„) = 0, Er, ^ fi'ki/2m 



j = 0,±l,±2,. 



(3.12) 



a well known result, which can also be obtained by solving the Schrodinger equation (IT) by means of separation 
of variablesQ. The present derivation of this result demonstrates the equivalence of the BIE (3.6) and the stationary 



Schrodinger equation. Note that because J-j{k) = {—ly Jj{k) [formula 9.1.5 in Ref. |2^ all the roots corresponding 
to j 7^ are doubly degenerate. 

Method II. Rather than approximating the BIE in Fourier space one can approximate it in coordinate space, 
i.e., one can solve (3.3) and not (B^). For this purpose one proceeds in two steps. First, one approximates the 
boundary F by a polygon with TV vertices situated on F, as shown in Fig. 0. Denoting the segment between vertices i 
and i + 1 by Fi one can write F fv U^;^Fi, and the BIE can be approximated by a sum of integrals along the N sides 
of the polygon. In a second step, one replaces along each segment F^ the function u{s) = d^ipnif) by a constant Mj. 
The BIE is then replaced by 



N 
i=l 



: f ds{r)G{r,r'; 



k) 



(3.13) 



Equation (3.13) still contains the continuous variable r' which should be eliminated. For this purpose, let us 
denote the position vector of the vertex i (see Fig.||) by Si and the position vector of the middle point of F^ by 
fi = (si + Si+i)J2rr, Then, setting in ( 3.13| ) r' — Vj, j = 1,2, ... ,N, one arrives at the so-called Boundary Element 
Equation (BEEpM 



N 



^u, As,/ ^d£_G{r,+£,As,,ry,k) = 0, 



(3.14) 



where Asi = s^+i — s^. The above equation represents a homogeneous system of N hnear equations and can be 
written 



N 






The elements of the matrix Bij(k), up to an irrelevant constant factor, are given by [cf. Eq.(3.1)] 



B,j (k) = Asj / d^ H^q'> {k \r, -r,+ CAsj\) 



(1) 



(3.15) 



(3.16) 



In analogy to our previous approach, the (approximate) energy eigenvalues can be obtained (in terms of fc) from 

det[Bij{k)] = 0, (3.17) 

i.e., as the real roots of this equation. 

3 _ _ 2 




i+1 N 

Fig. 2. The billiard boundary F is approximated by a polygon with TV vertices. 



The matrix elements Bij in (3.16) are expressed as single integrals in contrast to the matrix elements Aij defined 



in (3^) which are expressed in terms of double integrals. As a result, Method II is computationally less demanding 
than Method I, but has nevertheless two unfortunate features. First, the evaluation of the diagonal matrix elements 
Bii requires special integration technique due to the (integrable) singularity of the Green's function at ^ = 0. Second, 
in contrast to Method I where the truncation of the exact, infinite matrix Aij (defined in the Fourier space) provides 
us with a natural cutoff wave vector K^, in case of Method II the relationship between a similar K^ and the degree 
of discretization of the boundary (in real space) is less obvious. 

It should be emphasized itbat truncation in Fourier space is not quite equivalent to truncation (discretization of 
the boundary) in real spacellj. As an empirical rule, if one wishes to calculate energy eigenvalues corresponding to 
k ^ Kc accurately, one must take at least a few (about ten) discretization points for each section of the boundary of 
length equal to the corresponding de Broglie wave length A = 2t:/ Kc- Thus, the number of discretization points N 
scales with the length £ of the billiard boundary and the wave vector Kc as follows 



£. 10 

TV ^ 10- = — (KcC) ~ KcC. 
A Zn 



(3.18) 



Accordingly, accurate calculations of energy eigenvalues corresponding to sufficiently large k values require a large 
number of discretization points N along the boundary F, a condition which leads to large matrices Bij and, since 
these matrices are dense, to undesirable computational efforts. 

Method III. The most widely used method for the jeyaluation of billiard spectra is based on a non-singular 
version of the BIE (|3.3[). A simple, but not entirely rigorousO, derivation of this method applies the normal derivative 



oper ator d^' — i^ir') ■ Vr' to both sides of Eq.( pTTO ) which, according to definition (3.5) and with boundary condition 
( |3.2|) leads to 



h^ r 



(3.19) 



The integral kernel on the RHS, indeed, is non-singular at r = r'. BIE ( 3.19| ) is a homogeneous integral equation 
with unknown u(r); the energy eigenvalues E are given by the zeros of the corresponding Fredholni determinant [cf. 
Eq.(|3.4|)], i.e., as the solutions of 



det 



S{r-r') + —d^^G{r,r';E) 



= 



(3.20) 



Taking into account the explicit form (3J) of the free particle Green's function, Eq. (|3.19| ) can be written [cf. 
Eq.(PD] 



; (r') = —— (p ds{r) cos (j){r',r) H[ '{k\r' — r\)u{r) , 



(3.21) 



where 



cos 0(r',r) = i>'(r') 



(3.22) 



is the cosine of the angle between the exterior normal vector to T at r' and the unit vector corresponding to r' — r. 
Note that for r = r' the abov e sca lar product vanishes and, actually, cancels the singularity due to the Hankel function 
in the integrand of the BIE ( |3.2l| ). 

For a billiard with arbitrary boundary, the above functional determinant cannot be calculated analytically and one 
needs to resort to a numerical solution. For this purpose, one e mploys the same stra tegy as in case of Method II. 
After discretizing the boundary F, one can replace the BIE (3.21) by the BEE [cf. Eq.(3.14)] 



ik 



N 



— ^UiAsi d^ cos<l){rj,r^+£,Asi)H[^'{k\rj -r^-£,Asi\) , 



(3.23) 



where the notations are the same as in the case of Method II. Since the integrands on the RHS of the above equation 
are well behaved for all i,j — 1, 2, . . . , iV, one can approximate the corresponding integrals by the trapezoidal rule. 
As a result one obtains the system of linear equations 



N 



Y,G,,{k)u, = 0, 
i=i 



where 



dj yk) 



ik 



Wi 



Asj cos 4)ij H\ {kr. 



(3.24) 



(3.25) 



r, - r,- 



(3.26) 



The (approximate) energy eigenvalues can be determined as the roots of the determinant of Cij 

det [C,j{k)] = 0. 



(3.27) 




Fig. 3. Plot of abs_det(k) (open circles) for the circle billiard, corresponding to 60 and 300 discretization points (dp) of the 
boundary, by using Method III. The positions of the minima approximate the sought eigenvalues fc„. While the values of the 
minima depend strongly on the degree of discretization of the boundary, the actual positions of the minima do not. 

TABLE I. The first 18 distinct eigenvalues fc„ corresponding to the circle billiard of unit radius obtained by using Methods 
I, II and III. The eigenvalues corresponding to Method I represent the zeros of the integer Bessel functions Jj{k) and should 
be regarded as exact solutions. In case of Method II (III) the boundary was discretized by employing 60 (300) equally spaced 
points along the circle. The relative error for each approximate solution is less than 0.1%. 
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10.17347 
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7.0160 


11.08637 


11.0945 


11.0865 


7.58834 


7.5960 


7.5888 


11.61984 


11.6335 


11.6200 


8.41724 


8.4265 


8.4175 


11.79153 


11.8055 


11.7920 


8.65372 


8.6640 


8.6545 


12.22509 


12.2320 


12.2265 



B. Numerical Algorithm for Solving the BEE 



Based on the computational methods introduced we have written a FORTRAN _Z.7 program which implements the 
necessary algorithmic steps using the SLATEC Common Mathematical Library^. For all three methods one can 
employ a common algorithmic framework containing (i) a function det(fc) which, for an input wave vector fc, returns 
the complex value of the determinant of the corresponding system matrix, i.e., Aij, Bij or C^; (ii) a routine solve 
which calculates approximately the roots of the equation det(fe) = 0. Once the function det(fc) and the corresponding 
root finder solve are available one can scan the interval of k values of interest (between zero and the cut-off wave 
vector K^) to determine the zeros fc„ of det(fc) and, hence, the energy eigenvalues i?„. The algorithm may fail in 
practice when the separation between two consecutive eigenvalues is smaller than the scanning step Afc, i.e., when 
eigenvalues are nearly degenerate. The only way to avoid this error is to use the smallest affordable A/c. 

The actual form of det(fc) depends on the method chosen. In case of Method I, each matrix element Aij is given 
by a two-dimensional integral [see Eq.(3.9)] with singular and oscillatory integrand such that the evaluation of det(fc) 
would be computationally extremely demanding and would require special integration routines. Hence, we did not 
pursue an implementation of det(fc) for Method I. In the case of Methods II and III the function det(fc) consists of 
the following three parts 



(i) The subroutine discretize which takes as input the data necessary to define the actual form of the biUiard 
boundary and the number of discretization points N of the biUiard boundary; discretize returns as output 
the vectors r^, Si (i — 1,.. . ,N) [see Fig. |[ and other useful quantities based on them, such as the matrix 
'I'ij = fi — Tj, the vectors As^ = s^+i — s^, As^ = |Asi|, i>i = z x Asi/Asi (i.e., the external unit vector to the 
boundary at r^), etc. If one does not want to change the degree of discretization of the billiard boundary during 
the successive evaluations of det(fc), subroutine discretize should be run only once, namely during the first 
call of the function det(fc). 

(ii) The subroutine sysjnat which evaluates the complex valued matrix elements Bij and C'ij in-case of Method II 
and III, respectively. The Bij are evaluated according to Eq.(3.16) employing two SLATECEj (more precisely 



QUADPACl£3) quadrature routines, namely DQAGS, for calculating the non-diagonal matrix elements, and 
DQAWS, for calculating the diagonal matrix ele ments in w hich the integrand has a logarithmic singularity at 
^ = 0. The Cij are evaluated according to Eqs.( 3.2q - 3.26 ) in a straightforward way. In both cases the Hankel 



functions can be expressed in terms of the corresponding Bessel and Neumann functions for which the double 
precision SLATEC routines DBESJO, DBESJl and DBESYO, DBESYl, respectively, are called. 

(iii) The function det(fc) which evakiatcs the determinant of Bij and Cij, respectively. For this purpose one em- 
ploys the SLATEC subroutinesEj ZGEFA (factors a complex matrix by using Gaussian elimination) and ZGEDI 
(calculates the determinant and the inverse of a complex matrix by using the factors from ZGEFA) . 

The function det(fc) is complex- valued and, therefore, its real roots fc„ (the sought eigenvalues) must be simulta- 
neously zeros of both real and imaginary parts of this function. Due to the finite discretization of the boundary, the 
numerical solutions of the equation det(fc) = will be complex with a (hopefully) small imaginary part. In fact, 
the magnitude of the imaginary part of the "complex eigenvalue" fc„ can be used to characterize the accuracy of the 
energy eigenvalues thus determined through the real part of fc„. To the best of our knowledge, there exists no public 
domain subroutine which calculates automatically the roots of an arbitrary complex function of one complex variable, 
and as a result one can make little or no progress at all in the endeavor of constructing a robust fc„ eigenvalue finder 
algorithm based on the above straightforward approach. Howp^ter, there is a relatively simple solution to this problem 
which seems to be widely used by practitioners of the BIMIlilO. One notes that the zeros of det(fc) are also absolute 
minima for the square of the absolute value of this function, i.e., of abs_det(fc) = Re[det(fc)]^ -I- Im[det(fc)]^. Strictly 
speaking, the minima should assume zero values. The discretization of the boundary (or equivalently, the truncation 
of the original functional determinant) introduces errors such that the numerically evaluated minima of abs_det(fc), 
are small, but not zero; the magnitude of each minimum can be used to distinguish a real root of det(/c) from a local 
minimum of abs_det(/c). Since, numerically, it is much easier to determine the (local) minima of a real function of a 
real variable than to determine the roots of a complex function of a complex variable the suggested approach is much 
more convenient for our purpose. Accordingly, solve determines actually the local minima of abs_det(fc) by going 
through a given interval of wave vectors kmin l£ k < kmax in steps of Afc. Once a minimum is bracketed, its actual 
value can be calculated with any desired accuracy (for a given degree of discretization of the billiard boundary) by 
employing, for example, a double precision version of the function brent from Ref. ^ 

C. Numerical Results 



As a test of the algorithms described in Sec. |IIIB| and their numerical implementation w e det ermine the spectrum 
of a circular billiard. In this case Method I yields the exact energy eigenvalues fc„ [cf. Eq. ( 3.12| )] as the roots of the 



integer Bessel functions (these roots are in fact tabulated; see, e.g., Ref. 23). The first 18 distinct eigenvalues fc„ were 



also determined by means of Methods II and III described in Sec. [II A and are compared in Table | with the results of 
Method I. In case of Method II (III) 60 (300) equally spaced discretization points of the circular boundary have been 
employed. The locations of the minima of the function abs_det(k) have been determined by scanning the 2 < fc < 13 
interval with a step Ak = 0.004. Figure |3| illustrates the fc-dependence of abs_det(k), evaluated in the framework of 
Method III for two different discretizations of the boundary. An increase of the number of discretization points from 
60 to 300 changes significantly the values, but not the positions of the minima and, hence, does not affect significantly 
the values kn. 

Table || demonstrates that the res ults o f both Methods II and III reproduce the exact eigenvalues to at least three 



significant digits for k ^ Kc [cf. Eq. ( 3.1S|)] . In case of Method III, we have found that 300 discretization points lead 



to a precision of better than 1% for the 150 lowest eigenvalues of the circular billiard (with unit radius) corresponding 
to fc < 35. For larger k values the density of eigenvalues increases and, in order to separate adjacent minima of 
abs_det(k), one needs to reduce the step size Afc. The method breaks down for fc ^ C/N , i.e., when the distance 

10 



even-even 



odd-even 



even-odd 



odd-odd 



Fig. 4. The stadium billiard is formed by two semicircles of radius R connected through two parallel linear segments of length 
L. The two symmetry axes of the stadium are labeled H (horizontal) and V (vertical), and the corresponding four symmetry 
classes of the energy eigenfunctions are shown. 




between two consecutive discretization points of the boundary becomes comparable with the de Broglie wavelength 
of the particle, and the only remedy is to increase N . 
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Fig. 5. Spectral staircase A'^(fc) for the lowest 50 (70) energy levels of the (a) circle and (b) stadium billiard. In (a) the 
dashed line corresponds to the leading semiclassical Weyl term {Sk^ /A-k)/2, where the extra factor of 1/2 accounts for the 
double degeneracy of the energy eigenvalues with m 7^ 0. The solid line is obtained by taking into account the perimeter (next 
to the leading) term in the Weyl formula which for the circle billiard is given by AA^ = fc/4. In (b) the solid line corresponds 
to the asymptotic Weyl formula with the perimeter correction term. 



The program implementing Methods II, III allows one to calculate the spectra of billiards of arbitrary shapes, for 
which purpose one needs to solely alter the coordinates of the discretization points of the billiard boundary. As an 
example, we choose the Bunimovich stadium billiard depicted in Fig. which consist of two semi-circles (of radius 
R = 1) connected by two parallel linear segments (of length L) . We seek to calculate the lowest few hundred energy 
eigenvalues of both the circle and the stadium billiard. 

The circle billiard constitutes an integrable system, i.e., the number of constants of motion (energy and angular 
momentum) is equal to the number of degrees of freedom . Its energy eigenstates can be classified according to 
symmetry, i.e., by an orbital quantum number m, which counts the nodal lines through the center, and a principal 
quantum number n, which counts the nodes of the radial wave function, i.e., the nodal circlesQ. In contj|ia*t, the 
stadium billiard, regardless of how small L is, constitutes a non-integrable, i.e., (strongly) chaotic, systemErEZI. The 
study of quantum systems for which the underlying classical motion is chaotic is a relatively new and still widely open 
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field of studyH§E3. Since it is beyond the scope of the present article to provide an introduction to quantum chaos, 
we will content ourselves with considering without explanation one characteristic which distinguishes the spectra of 
non-chaotic (e.g., of a circle bilhard) and of chaotic (e.g., ofj-a jStadium biUiard) quantum systems, namely the so 
called (energy) level spacing distribution P{s). By definitionEj'u, P{s)ds is the probability that, given an energy 
level at E, .the nearest neighbor energy level is located in the interval ds about E + s. According to Random Matrix 
Theorip^T^ (RMT), applicable due to a quasi-random character of the Hamiltonian matrix, quantum systems, as 
far as the statistics of their energy spectrum is concerned, injgcncral can be classified into four universality classes, 
with well defined and distinct P{s) level spacing distribution£j'EJ. Integrable systems are described by the Poisson 
distribution with 



Pais) 



(3.28) 



The energy levels of classically chaotic systems, which do not break time reversal symmetry, (e.g., the stadium billiard) 
form a Gaussian Orthogonal Ensemble (GOE) with 



PcoEis) 



-sexp 



ns 



(3.29) 



Further universality classes are the Gaussian Unitary Ensemble (GUE) and the Gaussian Symplectic Ensemble (GSE); 
classical chaotic systems which break time reversal symmetry, e.g., ellipse or stadium billiards in an external magnetic 
field, belong to the GUE, while classical chaotic systems which preserve time reversal symmetry, but break spin 
rotational symmetry, e.g., a chaotic billiard in the presence of spin-orbit interaction, belong to the GSE. 



(a) 



(b) 
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— 4 independent GOE' s 







s s 

Fig. 6. Histogram of the energy level spacing distribution P{s) for the (a) circle and (b) stadium billiards. In (a) the histogram 
is constructed from the lowest 1,200 energy levels of the circle billiard. The solid line corresponds to the Poisson prediction 
for the level spacing distribution. In (b) the lowest 600 energy levels have been used to construct the histogram. The dotted 
(dashed) line represents the Poisson (GOE) prediction for P{s). The histogram is best approximated by the superposition of 
4 independent GOE distributions (solid line) which correspond to the same number of distinct symmetry classes of the energy 
eigenstates in a stadium billiard. 

Poisson and GOE distributions are distinguished most clearly near s = 0, since Po{0) = 1 constitutes the maximum 
of Po while Pgob(O) — constitutes the minimum of Pgoe', neighboring energy levels are likely to attract (repel) 
each other in the case of integrable (chaotic) systems. We want to show that the level spacing distribution evaluated 
by means of the BIM for circle and stadium billiards satisfies the Poisson and GOE distribution, respectively. For 
this purpose one needs to calculate at least a few hundred of the lowest energy levels without actually missing any 
energy levels since such misses would distort the energy level spacing distribution. The quality of the calculations, 
in particular in the case of the stadium billiard, can be judged from a comparison of the obtained (energy) staircase 
function M{E) (wiiici!. gives the number of quantum states with energy less or equal to E) with the corresponding 
Weyl-type formulaEila 
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Fig. 7. Histogram of the energy level spacing distribution P{s) for the (a) quarter- and (b) half-stadium billiards. The level 
spacing distribution for a quarter- (half-) stadium is well approximated by the GOE (two independent sequences of GOE) 
distribution function. 



(■^(^)> = J^iSE 



Cv^ + C 



(3.30) 



where S and C are the area and perimeter of the billiard, and C is a constant related to the geometry and topology 
of the billiard boundary. Presently, we employ units in which h /2m is equal to one; thus, e.g., E = k^ . Also, in the 
numerical results report ed be low we have chosen L ^ R ^ 1 (see Fig. m. 

Strictly speaking Eg. ( 3. 30 ) is valid only in the semi-classical (E — > oo) limit, but in practice it turns out that one 
can apply Weyl's formula even at the lower end of the energy spectrum. Our results for the staircase function A/'(fc), 
corresponding to the first 50 (70) distinct energy levels of the circle (stadium) billiard, are presented in Fig. ^. In the 
case of the circle billiard a complication arises due to the fact that all the energy levels with angular momentum m ^ 
are doubly degenerate. A simple remedy to this problem is to assume that the fraction of energy levels corresponding 
to 771 = is negligible in com parison to those with m ^ 0, and that the double degeneracy can be accounted for by 
dividing the RHS of Eq.(|33o|) by two. 

Based on the good agreement between J\f{E) and {Af{E)) shown in Fig. |^, one may conclude that all energy levels 
have been accounted for. A similar analysis for the first 600 energy levels showed that at most a few percent of the 
energy levels might been missed. This conclusion is independent of the method chosen, i.e., of Methods II and III. 

For a proper analysis of the energy level statistics we linearly scale the set of energy eigenvalues such that for the 
resulting sequeace^the mean level spacing is uniform and equal to unity. This transformation, known as "unfolding 
the spectrum"E3'til, is commonly achieved by replacing the original set of eigenenergies En = fe^ by 



Er, 



mEn)) 



(3.31) 



The unfolded spectrum now can be used to calculate the nearest level spacings s„ = En+i — En, which fluctuate 
about their mean value equal to one. Finally, a normalized histogram of the s„ series gives a rough representation 
of the distribution function P{s). The resulting distributions P{s) for the circle and stadium billiards are shown in 
Fig. p. In the case of the circle billiard the obtained histogram agrees very well with the expected Poisson distribution 
Eq.( |3.28 ). However, in the case of the stadium billiard the P{s) histogram does not resemble a GOE distribution, in 
particular, the distribution exhibits a clear absence of level repulsion, i.e., P{s) does not vanish for s — > 0. 

The deviation of P{s) from a GOE distribution arises due to the fact-that the stadium billiard, even though 
it is chaotic, exhibits a geometrical symmetry with two symmetry planed, as shown in Fig. ^ Accordingly, the 
stationary states fall into four distinct symmetry classes, according to their parity (i.e., either odd or even) with 
respect to reflection at the two planes. As a result, the stadium billiard spectrum is composed of four independent 
family of states, each of which is expected to conform to a GOE distribution. A general expression for the level spacing 
distribution function p(^'(s) corresponding to the superposition of N independent spectra with GOE statistics is 
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TABLE II. Wave vector eigenvalues fc„ (< 10) for (a) quarter-, (b) horizontal half-, (c) vertical half- and (d) full stadium 
billiards. The quarter stadium has only odd-odd energy eigenstates, the horizontal (vertical) half stadium has both odd-odd 
and odd-even (even-odd) eigenstates, while the full stadium has eigenstates which belong to all four symmetry classes. The 
symmetry of each eigenstate of the (full) stadium billiard can be identified based on this table as explained in the text. The 
dash in each column indicates that the corresponding eigenvalue is absent for that system. 



quarter 
stadium 


horizontal half 
stadium 


vertical half 
stadium 


full 
stadium 


quarter 
stadium 


horizontal half 
stadium 


vertical half 
stadium 


full 
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- 


- 


2.7784 


2.7785 


7.5231 


7.5238 


7.5238 


7.5240 


- 


3.4037 


- 


3.4037 


- 


7.6642 


- 


7.6640 


- 


- 


- 


3.7211 


- 


- 


- 


7.9760 


4.0564 


4.0565 


4.0565 


4.0566 


- 


- 


- 


8.0945 


- 


- 


4.6786 


4.6785 


- 


- 


8.3192 


8.3200 


- 


4.8800 


- 


4.8800 


- 


8.3989 


- 


8.3985 


- 


- 


- 


4.9223 


8.4639 


8.4642 


8.4639 


8.4640 


- 


- 


5.4931 


5.4935 


- 


- 


8.5200 


8.5200 


- 


- 


- 


5.6360 


- 


- 


9.0100 


9.0105 


5.7456 


5.7456 


5.7456 


5.7455 
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- 


- 


9.0600 


- 


- 


- 


6.2714 


9.2641 


9.2650 


9.2650 


9.2655 


- 


6.4387 


- 


6.4385 


- 


9.2890 


- 


9.2895 


- 


- 


6.5743 


6.5751 


- 


- 




9.3200 


- 


6.6493 


- 


6.6495 


- 


9.5900 


- 


9.5895 


6.9526 


6.9531 


6.9526 


6.9528 


- 


- 


9.8281 


9.8280 


- 


- 


7.1350 


7.1352 


9.9481 


9.9481 


9.9481 


9.9480 


- 


- 


- 


7.4815 


- 


- 


- 


9.9720 



derived in Appendix D. 
Eq.(pl) with AT = 4, i.e. 



Thus, the level spacing distribution corresponding to the full stadium billiard is given by 
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where erfc(z ) is t he complementary error functioned. Comparison of the numerically determined P{s) with the 
distribution ( 3.32|) in Fig. ^ is indeed satisfactory. The small values of P{s) for small s- values, i.e., values below the 
prediction by ( ^.32 ), are likely due to an omission of "nearly degenerate" eigenvalues by our spectrum finder routine 
(see also below). 

In order to check the assertion made about the symmetry classes of the energy eigenstates, and about the corre- 
sponding level spacing distributions, we have calculated and analyzed also the energy spectrum of a quarter stadium, 
and of the upper (horizontal) half and right (vertical) half stadium billiards, as well. The results are shown in Fig. [^. 
Indeed, the P{s) histogram for the quarter stadium, which accommodates all the eigenstates with odd-odd symmetry 
(see Fig. Q) conforms to a GOE distribution. On the other hand, for each of the two half stadiums, with eigenstates 
which belong to two distinct symmetry classes, namely odd-odd and odd-even (even-odd) in the case of horizontal 
(vertical) half stadiums, the level spacing distribution histogram is in good agreement with the theoretical prediction 
of the superposition of two independent GOE's as described by Eq.(D6) with N = 2, i.e.. 
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It should be noted that one can also identify the symmetry of each energy level of the stadium billiard. For this 
purpose one needs the energy spectrum of the full-, quarter-, horizontal half- and vertical half stadiums. These 
eigenenergies, corresponding to /c„ < 10, are listed in a convenient way in Table |lj The odd-odd eigenvalues can 
be simply read out from the column which contains the spectrum of the quarter billiard. Obviously, this eigenvalues 
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belongs also to the other three billiards under consideration. The odd-even (even-odd) eigenvalues can be obtained 
from the spectrum of the horizontal (vertical) half stadium by removing from the corresponding spectrum all the 
already known odd-odd eigenvalues. Finally, all the energy levels of the full stadium which have not been accounted 

for so far have even-even parity. 

We conclude this section with a few comments on the distribution function [Eq. (Pq)] 



P(^)(s) = 



92 



erfc 



/TT S 

2~iV 



N 



describing the superposition of N GOE distributions. For iV = 1 one recovers the GOE distribution function ( 3.29| ) 
wich is nownalized and yields a mean level spacing equal to onCj In the limit A'^ — > cxd, by using the the series 



expansiorO erfc(z) = 1 — {2/^/t:)z + O (a;'^) and the definitioiiE3 exp(— x) = limAr_»oo(l — x/N)^, one arrives at 
p(°°'>(^s) — exp (— s), which is exactly the Poisson distribution ( 3.28 ). This result is a particular case of the theorem 
according to which the level spacing distribution of the smwposition of infinitely many independent spectra (with 
arbitrary level spacing distributions) is always Poisson liktKI^Eil. 




Al' 



Fig. 8. (a) Geometrical construction of the deformed billiard (A'^ = 3), starting form the unit circle billiard, (b) The stadium 
billiard as an A^ = 2 deformed billiard, (c) Deformed billiard for N — 5. The highlighted regions correspond to "elementary 
sectors" for which simple GOE level spacing distribution is expected. 

Inspired by the billiard stadium problem, we propose a closely related numerical experiment which tests the appear- 
ance of the distribution P^-^\s) (D6). For this purpose we determine the energy levels corresponding to a deformation 
of the circle billiard involving an iV-fold symmetry axis. Let us consider N equidistant points Ai, i = 1, ...,N on 
the unit circle, with center O. Bi is the midpoint of the arc of circle AiAi^i. We construct then points Ci and Di 
by translating Bi with the vectors e • OAi_|_i and e • OAi, respectively. The parameter e controls the degree of the 
deformation. The deformed billiard is defined by the linear segments DiCi and the arcs of circle CiDi+i with unit 
radii. The new billiard, for A^ = 3 and e— 1 , is illustrated in Fig. ga. In the limit e — > one recovers the original circle 
billiard. For N — 2, the new billiard is actually the stadium billiard, as shown in Fig. g|b. For N > 2, the billiards 
can be regarded as a generalization of the stadium billiard; this is illustrated for another case, N ^ 5, in Fig. He. 

For a given A^, the deformed billiard possesses A^ symmetry planes and, therefore, the corresponding stationary 
states fall into 2A^ distinct symmetry classes, according to their parity with respect to reflection at these A^ planes. 
Proceeding as in the case of the stadium billiard, one can divide the deformed billiard into 2N elementary sectors 
(see the highlighted regions in Fig. S). The energy levels of a single sector should have an energy spectrum with 
GOE statistics. This is, indeed, born out of a BEM calculation as shown by the corresponding match with a GOE 
distribution in Fig. 0a in case of a single N — 5 sector. A billiard formed by attaching n {n < 2N) such elementary 
sectors should exhibit a level spacing distribution given by P^"''{s), while the level spacing distribution corresponding 
to the full deformed billiard should conform to P'-^^^{s). The level spacing distribution of an A^ = 5 billiard conforms 
well to the distribution P(^°)(s) as seen in Fig. 0b. 
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Fig. 9. Histogram of the energy level spacing distribution P{s) corresponding to A'^ = 5 for (a) an elementary sector, and (b) 
full deformed billiards. The level spacing distribution for the elementary sector (full) deformed billiard is well approximated by 
the GOE [-?(«)] distribution function (solid line). P{s) for the full billiard differs only slightly from the Poisson distribution 
(dashed line). 

In the limit iV — > cxd the deformed billiard becomes a circle of radius 1 + e as one can infer readily from the 
construction presented in Fig. 0. The suggested billiards produce then level spacing distributions which, due to 
limjv— ►oo P^^'{s) = Po{s), conform to a Poisson distribution. This is to be expected, of course, since this distribution 
governs the spectrum of a circle billiard. One can recognize in Fig. 0b that already in the case N ^ 5 the level 
spacing distribution resembles the Poisson distribution. 

Many further billiards can be constructed in a similar way. For the case of classical systems, a family of billiards 
which exhibit chaotic as well as mixed chaotic and regular motion have been studied in Ref. p5l The application of 
the BEM to determine level statistics as well as wave functions for the mixed system might reveal some surprising 
behavior. 

IV. OTHER EXAMPLES 

In this section we wish to present two other examples in which the BIM can be applied. Both examples exhibit 
the features mentioned at the end of Sec. O: (i) the corresponding Green's function is known analytically; (ii) the 
boundary condition at F assumes a simple form. Due to lack of space we shall only outline the BIM treatment of 
these examples. The interested reader is encouraged to work out further details, including the statistical analysis of 
the obtained data, in analogy to the quantum billiard case presented in the previous section. 

A. Finite Potential Well 

As a first example let us consider a particle trapped inside a two-dimensional potential well defined by a finite 
potential increase at the boundary, described by the potential 



T/(r) 





Vo 



for r G Vi 
for r E Do 



(4.1) 



Here 25^/0 represents the inner/outer region determined by a closed boundary F of arbitrary shape. The depth of the 
potential well is Vo {> 0). The energy spectrum for this system has a discrete part for En < Vo, and a continuous part 
for E > Vo- The quantum billiard studied in the previous section can be regarded as a limiting case of the present 
case corresponding to Vo ^ oo. 
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For the purpose of calculating the discrete energy eigenvalues of the system one applies the BIM presented in Sec. || 
for both inner (Vi) and outer (I?o) regions. As a result one obtains a set of two coupled BIE's; the two unknown 
functions are the wave function ^n and its outward (with respect to the inner region Vi) normal derivative di,ijjn 
along r. 



For V = 'Di, in analogy to Eq.(2.1C), the corresponding BIE reads 



i'^^Hr') = —P ids{r) UW(r)a,G« (r,r';i?„) - G« {r , r' ; E^,) d,^^\r) 
m Jr ^ 



(4.2a) 



with the Green's function 



G«(r,r';ii;„) = --^i/W (fc„ |r- - r' 



2r 



y/2mEn/h . 



(4.2b) 



The "exterior problem" T) = Do requires a more careful treatment due to the fact that I?o is unbounded. One can 
circumvent this difficulty by considering instead a finite region Vp delimited by F inside and by a circle Cp with a very 
large radius p outside, the center of the later located somewhere inside the regi on V j; evidently, limp^oo 2?p = T^o- 
Thus, when we apply Green's formula to obtain the BIE an extra term results in ( 2.1C ) due to the circle Cp. However, 
this extra term vanishes in the limit p ^ oo because for bound states (the only ones we are interested in) both the 
wave function i/jn and its gradient VV'n vanish exponentially at infinity. Hence, the corresponding BIE becomes 



^k°\r') = -^vids{r) Ui°Hr)a,G(°)(r,r';i?„)-G(°)(r,r';ii;„)a,Vi°'(0 

with the Green's function (which is finite for \r ~ r'\ -^ cx)) 

G(°)(r,r';i?„) = -l^ijW (»g„ |r - r'|) = --Ko{qn\r - r'\) , 
2a tt 



(4.3a) 



(4.3b) 



g„ = ^/2m{Vo-Er,)/h . 

Her e Ko {z) is a Bessel function of imaginary argumentE3 (see also Appendix A). Note the minus sign on the RHS of 
Eq.(4.3a) which accounts for the opposite orientation of the exterior normal unit vectors corresponding to 2?j and Vq- 
Since the wave function and its normal derivative must be continuous across F, i.e.. 



do) 



r^'ir) = rn">{r) ^ Mr) 



(4.4a) 



d.^i:\r) = d.i^l^Hr) ^ d.Mr), 



(4.4b) 



Eqs.(4.2-4.3) form a set of coupled BIE's with respect to the unknown functions ^„ and d^ipn- 

The numerical calculation of the energy levels of a particle in a finite two-dimensional potential we ll proceeds 
similarly as in the case of a quantum billiard. The steps to be filled in are the same as those discussed in Sees. Ill A HI B . 
Note, however, that due to the simultaneous presence of both ipn and du^Jn in the BIE's, only Method II can be applied 
in this particular case. 



B. Quantum Billiard in a Magnetic Field 



As a second example, let us consider a charged particle confined to a two-dimensional billiard with Vo — >■ oo [cf. 
Eq.(4.1)], in the presence of a uniform magnetic field B per pend icular to the plane of motion. The Hamiltonian for 
such quantum billiard in a magnetic field is given by [cf. Eq. (^.l|)] 



2m 



(4.5) 



where p = —ihW is the momentum operator, q is the electric c harg e of the particle, A{r) is the vector potential 
{B — W X A) and V{r) is the scalar potential as given by Eq.(O). The energy spectrum of the system can be 
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determined by solving the Schrodinger equation (|l.l|) subject to the Dirichlet boundary condition (3.2). To derive 
the corresponding BIE we rewrite the Hamiltonian ( [4.5|) , recaUing that for a static magnetic field V ■ A — 0, 



H = V' 

2m 



2m m 



and define the Green's function G{r,r';E) as the solution of 

E- H*{r)\G{r,r':E) = S {r - r') , 



(4.6) 



(4.7) 



where H* is the complex conjugate_of the Hamiltonian (4.6). Note that H* =/= H, which implies that the magnetic 
field breaks time reversal symmetryEj. 

Using the same strategy as in Sec. O, one can derive the following BIE 



4-2 p 

Mr') = —vids{r)[Mr)d,Gir,r';En)-G{r,r';En)d,^n{r)] 
m Jy- 

+ i'^lds{r)A,{r)G{r,r'-E^)Mr) , 
m Jr 



(4.8) 



where A^, = v ■ A. Since the wave function ?/;„ vanishes along the boundary of the billiard F [cf. Eq.(^^)] the last term 
in Eq.( |4.S| ) can be dropp ed. As a result, we obtain formally the same BIE as in the field-free case, namely Eq.( ^.3[ ), 
or equivalently Eq.( 3.2C| ). Hence, the energy levels of a quantum billiard in a magnetic field can be determined as 
described in Sec. III. The only difference is that, instea|d,of the free particle Green's function, the Green's function 
of a charged particle in magnetic field needs to be usedE3. Here we assume a vector potential corresponding to the 
symmetric gauge (i.e., A — Bx r/2) 



G{ry-E^) = -^(^'^"^'^)/''(-^)rQ-^)e-^/^C/Q-^,l,-) 



(4.9) 



(. — ^Jh/rau) is the so-calledj-Hiagnetic length, lo = qB/m is the cyclotron frequency, e = E/hui, ^= (r — r')^/2£^, 
r(a;) is the Gamma functioned and U{a,b,x) is the logarithmic confluent hypergeometric functioiiE3. The derivation 
of Eq. ([4.9| ) is beyond the scope of this article; the reader is referred to Ref. |3^. r— ■ 

The above Green's function can be evaluated numerically by employing the double precision SLATEC subroutincsEil 
DGAMMA, for the function T{x), and DCHU, for U{a,b,x). Since the evaluation of the Green's function and its normal 
derivative (which can be expressed analytically) is very time consuming in the presence of a magnetic field, it is 
recommended to apply Method HI for determining the energy spectrum of the system. 



V. CONCLUSION 



In this article we have attempted to provide a self-contained, tutorial like introduction to the Boundary Integral 
Method for calculating single particle energy spectra in two-dimensional nano-devices. The BIM is suitable whenever 
(i) a Green's function G is available analytically and (ii) the boundary condition at the boundary of the device is 
fairly simple. The method applies to arbitrary shapes of the boundary. 

As we have shown, the BIM can be successfully applied to calculate the energy spectrum of quantum billiards, 
allowing one to investigate the quantum signatures of chaos in these systems. The numerical accuracy of the BIM 
strongly depends on the degree of discretization of the billiard boundary. Unfortunately, by increasing the number of 
discretization points along the billiard boundary, the needed computational resources seem to increase more rapidly 
than the accuracy of the calculated energy levels. Since the number of the n eeded discretization points along the 
billiard boundary scales linearly with the cutoff wave vector K^, [see Eq.( 3.18 )], one can conclude that, in fact, the 
BIM allows one to calculate the lowest few hundred energy levels of any quantum billiard. The determination of 
higher energy levels, in general, becomes computationally too expensive. Needles to say, the other existing numerical 
methods for solving the Schrodinger equation present similar or even more stringent limitations and altogether they 
perform worse than the BIM. 

In conclusion, we would like to mention a few experimental confirmations of the energy level fluctuations of quantum 
billiards described in this article. The revived interest in studying billiard spectra, in the context of quantum chaos, 
has resulted in beautiful microwave experimentg^U'ES designed to test the theoretical predictions, based mainly on 
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random matrix theories. These experiments exploit the analogy between the Schrodinger wave equation of a quantum 
particle in an infinite two-dimensional potential well and the Helmholtz equation of the electromagnetic field in a 
resonant cavity. Thus, by microwave measurements in the range of 0-25 GHz frequency in quasi two-dimensional 
cavities shaped, e.g., in the form of a quarter stadium billiard, up to few thousands eigenfrequencies were meas 



in Refs. 37, fiSl and found in agreement with spectra obtained by employing the BIM. Microwave measurement^ 
resulted also in the direct observation of the eigenfunctions in microwave cavities of different shapes; the eigenfunctions 
were also fouud to be in agreement with descriptions by means of the BIM. A very recent microwave ("photon") billiard 
measuremento allowed for the first time the direct experimental study of the energy level statistics in the presence 
of broken time reversal symmetry; the level spacing distribution was found to conform to a GUE form. 
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APPENDIX A: FREE PARTICLE GREEN'S FUNCTION IN 2D 

In this appendix we derive the expression of the free particle Green's function G (r, r'; E) in two spatial dimensions. 
The corresponding expression in d-dimensions can be obtained in a similar fashion. 

The free particle Green's function is defined as the solution of the equation [cf. Eq. ( |2.2| )] 

E+^Vl^G{ry-E) = J(r-r') , 
{Vl + e) a (r, r'; k) ^ -^Hr - r') , (Al) 



where k = \J2E /rah is the wave vector of the particle of energy E^ and we have replaced the energy variable in the 
Green's function with fc, i.e., G(fc) = G{E). By changing variables R — r — r' , which is equivalent to moving the 
origin of the coordinate system to the point r', the above equation becomes 

{Vl + k')G{R;k) - '^5{R), (A2) 

n 

where G{R; k) = G{R, 0; k). The fact that Eq. (|A2|) does not contain r' and depends only on R is the consequence of 
translational symmetry. 



One can solve (A2) by of Fourier transform. Inserting the Fourier representations 
and 



G{R\k) = I -^-^G{q;k)eTcp{iqR) , (A3) 



S{R) = I ^exp(iqR) (A4) 



in Eq.(A2), and identifying the Fourier coefficients on both sides of the resulting equation, one arrives at 

2m 



q^ + k^)G{q-k) = — . (A5) 



Inserting G from (A5) into (A3) results in 



2m /■ (Pq exp{iqR) 

'W J (2^)2 fc2_^2 



G{R;k) = -^ I jt;^^:^ — — ■ (A6) 
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The two-dimensional integral is evaluated by using polar coordinates q — {q, 9) as follows 



1*00 7 1 p'Ztt 

G(R;k) = — ^ / —: TT—- / d9cxp(iqRcos6) 

ttTi^ Jo fc - Q 27r Jq 



(A7) 



The second integral on the right hand side is identified as one of the integral representations of the 0-th order Bessel 
function Jo(gi?) [cf. formula 8.4111. in Ref. g6| and one obtains 



G{R;k) = -—2 



°° qMqR) 
q K 



(A8) 



The integral on the RHS of ( |AS| ) is ill defined due to the singularity of the integrand at g = ±fc. However, the 
integral can be regularized by adding to k an infinitely small, positive imaginary part, i.e., k —>■ k + ie. In this 
case k^ ^> {k + ie)"^ = — (e — ik)"^, and according to the formula 6.5324 of Ref. Oq the integral in (AS) is equal to 
Ko{{e — ik)R), where Kq is the MacDonald (modified Bessel) function, which is finite as i? ^ oo. After taking the 
limit £ — ^ 0^, one obtains then 



G{R;k) = ^Koi-ikR) 



(A9) 



Note that for e < the above integral would be divergent for R -^ cx3. However, as long as we are not concerned 
with the R ^ oo behavior of G{R; k), the infinitesimal e caii-4)e chosen either positive or negative. The choice £ > 
is equivalent to the so-called Sommerfeld radiation conditiorE^. 

Finally, by using the identity Ko{z) = {in/2)HQ (iz) [formula 8.4071 in Ref. 26|, where z is an arbitrary complex 



r(i) 



number and i?o i^ the 0-th order Hankel function of the first kind, one arrives at the expression of the free particle 



Green's function given by (3.1) 



APPENDIX B: EVALUATION OF THE SINGULAR INTEGRALS IN THE BOUNDARY INTEGRAL 

EQUATION 



In order to calculate the LHS of Eq.(2_^), consi der first the case when the potential energy V{r) is zero and, 
therefore, the relevant Green's function is given by (3.1). Fcp-£ = |r — r'| -^ one can replace the Hankel function in 
the above equation by its limiting form for small argument^ 



G{r,r';E„ 



nh' 



■ ln(fc£) , 







(Bl) 



Next, we parameterize the arc of circle Ce through the angle if (see Fig.|i|) formed by the tangent ATi to F at Ae F (of 
position vector r') and the vector e. The angle (p assumes values between zero and 0{r') — T1AT2, i-e., the exterior 
angle made by the two tangents to F at A. If the contour F is smooth then the tangents coincide and 0{r') — tt. The 
arc element along C^ is ds{r) — edip. Since both ipn and d^^pn are finite, in the limit £ — > the quantities can be 
replaced in (2^) by their values at r = r'; one obtains then 



lim-— / ds{r)G{r,r';En)d^i}nir) 
e^o 2m Jc^ 



-— liin[eHke)]dir')d,Mr') 







(B2) 



The integral containing d,yG in ( pT9| ) can be calculated in a similar fashion. According to Eqs.(2.7) and (3.1) one 
can write successively 



d,G{r,r';En) = i/(r) • V 
imk 



2r 



v{r) 



im 
2h 
r ~ r 



^H^'^{k\r-r'\ 



r(i) 



(B3) 



H','{k\r-r'\) , 



where we used dH^ '{z)/dz — —HI '(z) [cf. formula 9.1.30 in Ref. g^. On the arc of circle Cg the dot product in 
(B3) is equal to one (see Fig.|l|); taking into account the limiting form of H^ for small argumental3, one can write 
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lim 



ds{r)%l;n{r)d^G{r,r';En) = -— V'n(r') lim 



2m 

27r 



)(r') 



e(i(/5 



2r 



2i 

Trfce 



(B4) 



VA.(r') 



For a smooth boundary F, where 6'(r') = tt, Eqs.(B2-B4) provide the result given in (2.9) 



For a finite potential energy V{r)^ in general, there is no simple analytical expression for the Green's function G 
and the validity of Eq.(2.9) is question able . However, by assuming on physical grounds that V{r) is finite for all 
r £ V, one can realize that the result ( |2.9D holds in this case too. Indeed, when e is small the potential energy is 
almost constant in the vicinity of r' (point A in Fig.|l|) and, therefore, one can approximate the Green's function with 
the corresponding expression valid for a constant V = V{r'). The approximation becomes exact in the limit e — > 0. 
But G for a constant potential energy has essentially the same form as for a free particle [Eq^S.l)] and, therefore, it 
has the same type of logarithmic singularity at r = r'. Since the actual value of the integral ( |2.9| ) is determined solely 
by the type of this singularity of the Green's function one may conclude that the result derived in this appendix holds 
in general. 



APPENDIX C: ANALYTICAL SOLUTION OF THE BOUNDARY INTEGRAL EQUATION FOR A 

CIRCULAR BILLIARD 



In this appendix we solv e analytically the BIE (B.6) for a circular billiard of unit radius. For the unit circ le C — 2t: 
and, according to Eq.(3.7b), one finds Kj = j, with j — 0, ±1, . . .. By using the Fourier representation ( p. 7a ) for u{s), 
the BIE becomes 



oo p2tt 

y Uj / ds G {s , s' ; k) exp{ij s) = 

_ „ "'0 



j = -oo 



The expression (A6) of the free particle Green's function, in the present case, can be written 



2m f d^q exp{iq[r{s) - r{s')]} 



(CI) 



(G2) 



7r?i' 



f°° qdq 1 f"^^ 

■ / —t: TT — / d9exp\iqcos(s — 9)]exp\—iqcos(s' — 9)] 

Jo «: - Q 27r 7o 



where q and 9 are the polar coordinates of the 2D vector q. Inserting (G2) in (CI) one obtains (the irrelevant constant 
factor 2m /h^ can be dropped) 

°° 1'°° ndn 1 /■^'^ 

E -^X 1^2^ I dOcxj^Hqcosis'^e)] (G3) 

1 f^"" 

X — / ds exp[iq cos{s — 9)] exp{ij s) = 0. 
27rJo 

By taking into account the integral representation of the integer Bessel function [formula 8.4111. in Ref. Eq], the 
integral over s in ( |G3| ) can be evaluated exactly as follows 



I r2TT r -^ .27r 

— / ds exp[iq cos{s — 9)] exp{ij s) — tj- ds exp{i[q cos{s — 9) + j (s — 9)]} 

X exp{ij9) — f Jj{q)exp{ij9) . 



2tt 



(G4) 



Similarly, the integral with respect to 9 in (C3) gives 



1 f^"" 

— / d9ex-p[—iqcos{s' — 9)]exp{ij9) — {—iy Jj{q)exp{ijs') . 
27r Jo 



(G5) 



With the last two results, Eq.(C3) becomes 
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y Uj exp(jjs') / dq 



k^ — q^ 



= 0. 



(C6) 



By employing formulas 6.535, 8.4061 and 8.4071 in Ref. Eq, the integral with respect to q can also be calculated 
exactly as follows 






dq 



q[JAi)r 



9' + {ikf 



= -I,{-ik)K,{-ik) = - J^{k)H]'>(k) 



tWi 



(C7) 



Here Ij and Kj are imaginary argument Besscl functions. 

Finally, the BIE for the circle billiard of unit radius can be written as 



^ UjeMi]s')Jj{k)Hf\k) = 



(C8) 



Since exp(ijs), j = 0,±1..., form a complete orthonormal set, the above equation tells us that the expansion 
coefficients should be all equal to zero. The eigenenergies correspond to those k values for which non trivial Wj's 
exis t. Re calling that the Hankel functions have no real roots, one obtains the same eigenenergy equation (3.12) as in 



Sec.IIIA 



By taking the Fourier transform of Eq.(C8) with respect to s', one can see that the matrix Aij{k) defined by (3.9) 
is indeed diagonal and 



(1)/ 



Ay(fc) cx Jj{k)HY>{k)5, 



(C9) 



APPENDIX D: SUPERPOSITION OF iV INDEPENDENT SPECTRA WITH GOE LEVEL SPACING 

DISTRIBUTION 

In this Appendix we derive the level spacing distribution function P^^'{s) corresponding tp-|the superposition of N 
independent spectra with GOE level spacing distribution. In general, P{s) can be expressedE3 



P{s) 






(Dl) 



where £{s), the so-called gap probability, gives the probability that the energy interval {E, E + s) lacks energy levels. 
Let us consider N independent (i.e., uncorrelated) sets of energy levels with GOE level spacing distribution 






TT / S \ ' 

4 \n) 



i = 1,, 



The probability density (p3) is normalized as follows 

/ dsP,{s) = 1, (s,) = / dssP,{s) 

Jo Jo 



.N 



N 



(D2) 



(D3) 



Note that the choice (s^) = N for each set of levels leads to a unit mean level spacing (s*^^-') for the spectrum 
comprising all N ene rgy s pect ra. 

According to Eqs.(Dl)-(D2), the individual gap probabilities can be expressed as 



1 /"OO /'OO 



(D4) 



2N 



dt exp (— t^ 



erfc 



2N 



where erfc(z) is the complementary error functioncJ. Since the energy spectra are uncorrelated, the gap probability 
of the combined spectrum is given by 
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N 



£^^\s) = l[£M 



erfc 



2 N 



N 



(D5) 



and, according to (Dl), the desired level spacing distribution function becomes 

N 



p'-Hs) ^ ^ 



erfcl— - 



(D6) 



Note that the above method of calculating P{s) is rather general and applies also when the independent spectra 
have arbitrary statistics. For example, one could calculate the level spacing distribution of the superposition of an 
arbitrary number of spectra, some of them obeying Poisson statistics and the rest GOE statistics. For further details 
the reader is referred to Refs. ^ |2|. 
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